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H ■ Abstract 

^ , We introduce the smt toolbox for Matlab. It implements optimized 

I storage and fast arithmetics for circulant and Toeplitz matrices, and 

^T) ' is intended to be transparent to the user and easily extensible. It also 

. provides a set of test matrices, computation of circulant precondition- 

^ I ers, and two fast algorithms for Toeplitz linear systems. 

^ 1 Introduction 

I arise naturally in a large number of applications, like medical imaging, remote 
^ I sensing, geophysical prospection, image deblurring, etc. Moreover, in many 
t:;-]- ■ real-world computations, the full exploitation of the structure of the problem 
0^ . is essential to be able to manage large dimensions and real time processing. 
^ ! In the last 25 years, a great effort has been made to study the properties of 
^ I algebraic structures and to develop algorithms capable of taking advantage of 
O I these structures in the solution of various matrix problems (solution of linear 
^ ■ systems, eigenvalues computation, etc.), as well as in matrix arithmetics, for 
^ ■ what concerns memory storage, speed of computation and stability. 
• ^ . In spite of the many important advances in this field, there is not much 
^ I software publicly available for structured matrices computation. On the con- 
ed I trary, most of the fast algorithms which have been proposed, and whose 
properties have been studied theoretically, exist only under the form of pub- 
lished papers or, in some occasion, unreleased (and often unoptimized) re- 
search code. This fact often force researchers to re-implement from scratch 
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algorithms and BLAS-like routines, even for the most classical classes of 
structured matrices. Anyway, this is possible only for those with enough 
knowledge in Mathematics and Computer Science, and totally rules out a 
large amount of potential users of structured algorithms. 

Matlab [T3] is a computational environment which is extremely diffused 
among both applied mathematicians and engineers, in academic as well as 
in industrial research. It makes matrix computation sufficiently easy and 
immediate, and provides the user with powerful scientific visualization tools. 
At the moment, besides the standard unstructured (or full) matrices, the 
only matrix structure natively supported in Matlab is sparsity. In sparse 
matrix storage only nonzero elements are kept in memory, together with 
their position inside the matrix. Moreover, all operations between sparse 
matrices are redefined to reduce execution time and memory consumption, 
while mixed computations return full or sparse arrays depending on the type 
of operations involved. 

Our idea is to extend Matlab with a computational framework devoted 
to structured matrices, with the aim of making it easy to use and Matlab- 
like, transparent for the user, highly optimized for what concerns storage and 
complexity, and easily extensible. We tried to follow closely the way Matlab 
treats sparse matrices, and for doing this, we used the Matlab object-oriented 
classes. Starting from version 5, in fact, it is possible to add new data types 
in Matlab (classes), to define methods for classes, i.e., functions to create 
and manipulate variables belonging to the new data type, and to overload 
(redefine) the arithmetic operators for each new class. 

At the moment, our toolbox supports two very common classes of struc- 
tured matrices, namely circulant and Toeplitz matrices. In writing the soft- 
ware our aim was not only to furnish storage support, full arithmetics and 
some additional methods for these structured matrices, but also to create a 
framework easily extendible, in terms of functions and new data types, and 
to specify a pattern for future developments of the package. So, a great effort 
was spent in the software engineering of the toolbox. 

Among the available Matlab software for structured matrices computa- 
tion, we mention the following Internet resources. Various fast and super- 
fast algorithms for structured matrices have been developed by the MaSe- 
team (Matrices having Structure) [21] , coordinated by Marc Van Barel at the 
Katholieke Universiteit of Leuven. A Toolbox for Structured Matrix Decom- 
positions ^26] has been included in the SLICOT package [21] , developed under 
the NICONET (Numerics in Control Network) European Project [25j. The 
Restore Tools [H] is an object oriented Matlab package for image restoration 
which has been developed by James Nagy and his group at Emory Univer- 
sity. The MOORe Tools [llj, an object oriented toolbox for the solution of 



2 



discrete ill-posed problems derived from |8j, provides some support for cer- 
tain classes of structured matrices, mainly Kronecker products, circulant and 
block-circulant matrices. 

Matlab implementations of various algorithms are also available in the 
personal home pages of many researchers working in this field. Many sub- 
routines written in general purpose languages, like C or Fortran, are also 
available. It is worth mentioning that there are plans to add support for 
structured matrices in LAPACK and ScaLAPACK; see P, Section 4]. 

The plan of this paper is the following. In Section [21 we describe in detail 
our toolbox, called smt (Structured Matrix Toolbox), its capabilities, the new 
data types added to Matlab and the functions for their treatment. Section [3] 
is devoted to some technical implementation issues, while in Section H] we 
describe possible future lines of development of this software package. 

2 The toolbox 

Once installed (see Section [3]), the toolbox resides in the directory tree 
sketched in Fig. [H The main directory contains a set of general purpose 
functions, described in detail in Section [231 the following four subdirec- 
tories: 

• Osmcirc and Ssmtoep, which contain the functions to create and ma- 
nipulate the objects of class smcirc and smtoep, i.e., circulant and 
Toeplitz matrices; 

• private, whose functions, discussed in Section [231 are accessible by 
the user only through the commands placed at the upper directory 
level; 

• demo, which hosts an interactive tutorial on the basic use of the tool- 
box. 

file system 
smt 



private @smcirc @smtoep demo 



private 

Figure 1: Directory tree of smt 
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Let us briefly explain how Matlab deals with new data types. When 
the user creates an object of class, say, obj, then the interpreter looks for 
the function with the same name in a directory called @obj, located in the 
search path. Similarly, when an expression involves a variable of class obj or 
a function is applied to it, the same directory is searched for an appropriate 
operator or function defined for objects of this class. 

Writing this software, we took great care in checking the validity of the 
input parameters, in particular for what concerns dimensions and data types, 
and in using an appropriate style for warnings and errors, in order to guar- 
antee the Matlab-like behaviour of the toolbox. As this requires a long chain 
of conditional tests, the resulting functions are often more complicated than 
expected (see for example the file mtimes.m in the directory @smcirc), but 
this does not seem to have a significant impact on execution time. 

Full documentation for every function of the toolbox is accessible via the 
Matlab help command and the code itself is extensively commented. Manual 
pages can be obtained by the usual Matlab means, i.e., 

help <func_nanie> for the functions in the main directory, 

help <class>/<func_name> where <class> is either Osmcirc or Osmtoep, 
help private/<func_name> for the functions in the private subdirectory. 

Notice that <func_name> may be Contents (except in conjunction with 
private), in which description of the entire directory content is dis- 

played. For example, the command 

help Osmtoep/Contents 

displays the list of all the functions, operators and methods for smtoep ob- 
jects (i.e., Toeplitz matrices), while 

help Osmtoep/mtimes 

gives information about the matrix product operator for Toeplitz matrices. 

We remark that the toolbox supports both real and complex structured 
matrices, since complex numbers are natively implemented in Matlab. It is 
also possible to manage sparse smcirc and smtoep objects. 

2.1 Circulant matrices 

A circulant matrix [5J of order n is a matrix C whose elements satisfy the 
relations 

^ij 1, . . . , fl, 

Ck—n Cfc, k 1,...,?7. 1, 
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e.g., for n = 4, 



C = 



Co 


ca 


C2 


Cl 


Cl 


Co 


Ca 


C2 


C2 


Cl 


Co 


C3 


C3 


C2 


Cl 


Co 



The main property of a circulant matrix is that it is diagonahzed by the 
normahzed Fourier matrix Tr,, defined by 

where r\ is any primitive complex n-th root of unity (i.e., r^^^ 7^ 1 for A; = 

27ri 

0, . . . , n — 1, and 77" = 1). We let 77 = u; := e~ and T = T^^. 
This allows us to factorize any circulant matrix in the form 

C = JTAJT*, (2.1) 

where 

A = diag(c'(l),C'(^),...,(7(u;"-^)), 

and 

c'(c) = Ec,r'= 

fc=0 

is the discrete Fourier transform of the first column of C . Given the definition 
of discrete Fourier transform adopted in the f f t command of Matlab, we have 

6 := diag(A) = fft(c), 

being c = (cq, Ci, . . . , c„_i)'' the first column of C . 

In smt, a variable C of class smcirc is a record composed by 4 fields. The 
field C.type is set to the string 'circulanf, and is a reminder, present in all 
smt data types, denoting the kind of the structured matrix. The first column 
of the circulant matrix C gives complete information about it, and is stored 
in C.c, while C.dim is the dimension n. The field C.ev contains the vector 
6 of the eigenvalues of C; it is computed when the matrix is created and 
updated every time it is modified. This means that the initial allocation of 
a circulant matrix, as well as some operations involving it, takes O(nlogn) 
fioating point operations {flops). 

For example, an object of class smcirc can be created specifying its first 
column, with the command 

C=smcirc([l;2;3;4]) 
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and it is visualized either as a matrix 



C = 



1 

2 
3 
4 



4 
1 
2 
3 



3 
4 
1 
2 



2 
3 
4 
1 



or showing its record structure 



C = 



smcirc object with fields: 
type: 'circulant' 
c: [4x1 double] 
dim: 4 
ev: [4x1 double] 

depending on how the configuration parameter display is set; see the func- 
tion smtconf ig in Section I2.4[ The structure of the object can also be 
inspected with the command get(C), independently on the configuration of 
the package. If the column vector passed to smcirc is of class sparse, this 
memory storage class will be preserved in C . c, but not in C . ev. 

All operations between circulant matrices have been implemented, when 
possible, by fast algorithms, meaning that they require a complexity smaller 
than the corresponding unstructured matrix operations. For example, when 
the user computes the sum of two circulant matrices with the command 
A=C+D, the function plus.m is automatically called, in order to sum the .c 
fields, and to update the .ev field of the resulting object, as follows 



To multiply a circulant matrix C times a vector x, we can exploit the fac- 
torization fl2.1l) to obtain 



where u o v = {uiVi, . . . , UnVn)^ denotes the Schur product of two vectors. 
This requires only 2 f f t's, since the vector 6 is stored in the . ev field of the 
corresponding object C. In a similar way, the first column of the product of 
two circulant matrices is evaluated by the inverse discrete Fourier transform 
of the product of the eigenvalues of the two factors. In all cases, the computa- 
tion is optimized in terms of complexity. After performing many operations 
which update the eigenvalues of an smcirc object, it may be advisable to 
recompute C . ev, to improve its accuracy; if required, this can be done by 

C=smcirc (C . c) , 



A.c = 



C.c + D.c 



A.ev = C.ev + D.ev. 




(2.2) 
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Operators and special characters 



plus 


A+B 


power 


k.~2 


uplus 


+A 


mldivide 


A\B 


minus 


A-B 


mrdivide 


A/B 


uminus 


-A 


Idivide 


A.\B 


mtimes 


A*B 


rdivide 


A./B 


times 


A.*B 


transpose 


A. ' 


mpower 


k~2 


ctranspose 


A' 



Table 1: Overloaded operators 



as the user is not allowed to directly modify the fields of an object belonging 
to an smt class. 

All the overloaded operators, or methods, for smcirc objects are coded in 
a set of functions, whose names (fixed by the Matlab syntax) are reported 
in Table [H, together with the equivalent Matlab notations. Each of these 
functions is called when at least one of the operands in an expression is 
of class smcirc; if the two operands are different smt objects, the method 
corresponding to the first one is called. The result is structured whenever 
this is possible. 

When an operation is performed between two circulant matrices, the com- 
plexity is not larger than 0{n\ogn) (for example in the matrix product), 
while it may be larger when one of the arguments is unstructured; e.g., the 
product between a circulant and a full matrix, which is computed by multi- 
plying the first operand times each of the columns of the second one, takes 
0{n^ \ogn) flops. So, if C and D are both smcirc objects and x is a vector, 

E=C*D produces a smcirc object, 

y=(C+D)*x returns a vector, 

F=C*rand(n) returns an unstructured matrix, 

and, in all cases, the fast algorithms implemented for the smcirc class are 
automatically used in the computation. 

The operations on smcirc objects which rely on the factorization (12.11) . 
and so exhibit a 0(n"logn) complexity, are the product, the power and the 
left/right division for matrices. Obviously, some operators do not involve 
fioating point computations at all, like the transposition or the unary minus. 

A trivial implementation of the algorithms is not sufficient to obtain a 
package which is both robust and transparent for the user. In fact, each 
function should be able to handle most of possible user's errors, and should 
replicate the typical behaviour of Matlab when any of the operands are scalars 
or empty arrays. As an example, we report in Algorithm [T] the structure of 
the plus.m function, which is called when the first structured argument in 
a sum is an smcirc object. As it can be seen, when the operands are of 
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Algorithm 1 (The plus.m function for smcirc objects) 

check vaUdity and dimensions of input arguments C and D 
deal with scalar or empty arguments 
if C is smcirc 

if D is scalar or smcirc 
the result is smcirc 
else if D is smtoep 

the result is smtoep 

else 

the result is full 
else {i'f C is not smcirc, then D is] 

if C is scalar 

the result is smcirc 

else 

the result is full 



different classes, the result belongs to the less structured class; e.g., circulant 
plus full is full, circulant plus Toeplitz is Toeplitz, etc. 

Many Matlab standard functions have been redefined for circulant ma- 
trices. They are listed and briefiy described in Table O among them, there 
are simple manipulation and conversion functions, like abs, double or full, 
some which return logical values (the isXXX functions), and a few which 
optimize some computations for smcirc objects, like det, eig or inv. The 
implementation of the last three functions is straightforward, as each smcirc 
object contains the eigenvalue of the circulant matrix in its .ev field. We 
remark that some functions require a larger complexity for a circulant than 
for a full matrix, like imag, because extracting the imaginary part of the 
entries of a circulant matrix requires to recompute its eigenvalues. 

The list in Table [2] is surely incomplete, since in principle all Matlab 
matrix functions could be overloaded for circulant matrices. We implemented 
those functions which we consider useful, leaving an extension of this list, if 
motivated by real need, to future versions of the package. It is sufficiently 
easy to add new methods to the class, since the user can start from an existing 
function, as a template, and then place the new file in the smt/Osmcirc 
directory. 

Let us add some comments on some of the functions listed in Table [21 
When adding a new class to Matlab, there are a number of functions which 
must be defined so that the class conforms to Matlab syntax rules. The 
get method allows to extract a field from an object, while display defines 
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Elementary math functions 



abs 


absolute value 


fix 


round towards zero 


angle 


phase angle 


floor 


round towards —oo 


conj 


complex conjugate 


ceil 


round towards +00 


imag 


imaginary part 


round 


round argument 


real 


real part 


sign 


signum function 


Basic array information 


size 


size of array 


get 


get object fields 


length 


length of array 


is empty 


true for empty array 


display 


display array 


isequal 


true for equal arrays 


Array operations and manipulation 


diag 


diagonals of a matrix 


reshape 


change size 


full 


convert to full matrix 


tril 


lower triangular part 


prod 


product of elements 


triu 


upper triangular part 


sum 


sum of elements 






Array utility functions 


double 


convert to double 


subs as gn 


subscripted assignment 


single 


convert to single 


subsindex 


subscript index 


isa 


true if object is in a class 


subsref 


subscripted reference 


isf loat 


true for floating point 


end 


last index 


isreal 


true for real array 








Matrices and numerical linear al 


gebra 


det 


determinant 


inv 


matrix inverse 


eig 


eigenvalues and eigenvectors 







Table 2: Overloaded functions 



how an object should be visualized on the screen; this can be customized 
in smt, as it will be shown in Section 12. 4[ Some other functions define 
the effect of subindexing on the new class. We let two of them, subsasgn 
and subsindex, just return an error code for an smcirc object, since we 
consider them useless for circulant matrices. The third one, subsref, is 
a function which allows to access a field (C.c) or an element (C(2,3)) of 
a circulant matrix, and to use typical Matlab subindexing expressions like 
C(:) or C(3:4, :). Notice that C(l:3,4:7) returns a Toeplitz matrix (i.e., 
an smtoep object; see Section [2^21) . while C( [1 ,3 , 5] ,6 : 8) returns a full 
matrix. 

The class smcirc includes two additional methods: smtvalid is a func- 
tion, called by other functions of the toolbox, which determines if an object 
is a valid operand in an expression, while smtoep converts an smcirc object 
into an smtoep one, as a circulant matrix is also a Toeplitz matrix. 
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2.2 Toeplitz matrices 



A Toeplitz matrix of order n is a matrix T whose elements are constant along 
diagonals, that is 

Tij i-ij 1, . . . , 

e.g., for n = 4, 

^0 t-i i-2 ^-3 
rp _ ^1 ^0 ^-1 

t2 ti to ^-1 

h ti to 

We introduced a class smtoep, for Toeplitz matrices, similar to the smcirc 
class. An smtoep object can be created by specifying its first column and 
row, for example with the command 

T=smtoep( [4:7] , [4:-l:l]), 

or giving only the first column, in which case the resulting matrix is Hermi- 
tian. Similarly to what happens to smcirc objects, an smtoep object can be 
displayed either as 

T = 

4 3 2 

5 4 3 

6 5 4 

7 6 5 

or 

T = 

smtoep object with fields: 
type: 'toeplitz' 

t: [7x1 double] 
diml: 4 
dim2: 4 
cev: [8x1 double] 

depending on the display configuration parameter; see smtconf ig in Sec- 
tion [231 

An smtoep object has two fields for the number of rows (T.diml) and 
columns (T . dim2) of the matrix, while T . t contains the data to reconstruct 
the matrix, namely the first row and column, in a form which is convenient 
for computation; in the above example, 

T.t = (l,2,...,7)^. 



(2.3) 



1 

2 
3 
4 
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The meaning of the T.cev field will be explained later in this section. 

The componentwise operators, like sum, subtraction, and the so-called 
dot-operators of Matlab, can be easily implemented for Toeplitz matrices, 
similarly to what has been done for the smcirc class. 

Regarding the matrix product, it is well known that a Toeplitz matrix 
can be embedded in a circulant matrix Ct'-, e.g., given the matrix (12.31) . we 
can write 





t-l 


t-2 


t-3 





t3 


t2 


tl 


ti 


to 


t-l 


t-2 


t-3 





t3 


t2 




tl 


to 


t-l 


t-2 


t-3 





ts 


h 


t2 


tl 


to 


t-l 


t-2 


t-3 








t3 


t2 


tl 


to 


t-l 


t-2 


t-3 


t-3 





t3 


t2 


tl 


to 


t-l 


t-2 


t-2 


t-3 





ts 


t2 


tl 


to 


t-l 


t-l 


t-2 


t-3 





ts 


t2 


tl 


to 



(2.4) 



So, to compute the product y = Tx by a fast algorithm, one can construct a 
vector X by padding x with zeros to reach the dimension of Ct, then compute 
Cyx by (12. 2p . and finally extract y from the first components of the result. 

The zero diagonals in (12.41) can be deleted, in which case the dimension 
of Ct is minimal, or "tight" : if T is m x n, then the "tight" dimension of Ct 
is m + T2, — 1. On the contrary, we can insert as many zero diagonals as we 
want. This may be useful, because the implementations of the f f t perform 
better when the length of the input vector is a power of 2. 

In smt both choices are available, and can be selected by editing the 
command smtconst; see Section 12. 4[ Although Matlab implementation of 
the f ft, namely FFTW [7], exhibits a very good performance also when the 
size of the input vector is a prime number, we observed that matrix product 
is generally faster if we extend the matrix Ct to the next power of 2 exceeding 
m + n — 1. 

A particular function has been created to speed-up Toeplitz matrix mul- 
tiplication. Thus, the command 

T=toeprem(T) 

pre-computes the eigenvalues of the matrix Ct, and stores them in the . cev 
field. This is done automatically when an smtoep object is allocated, and 
allows to perform only two f ft's for each matrix product, instead of three. 
The price to pay is that, like in the case of circulant matrices, some ele- 
mentary functions involving smtoep objects have a complexity larger than 
expected, as they need to compute the .cev field of the result. If this be- 
haviour is not convenient, the automatic call to toeprem can be disabled 
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by the smtconfig command (see Section [2 ■41) . and the user can either call 
toeprem when needed, or renounce to multiplication speedup. 

All the operators and functions of Tables [T] and [2] have been implemented 
for the class smtoep, with some differences. Unlike the circulant matrices, 
there is not a standard method to invert a Toeplitz matrix, or to compute 
its determinant or eigenvalues. On the contrary, various different algorithms 
are available and, probably, more will be developed in the future. For this 
reason the functions inv, det and eig, supplied with the toolbox, return an 
error for an smtoep object, and they are intended to be overwritten by user 
supplied programs. 

2.3 Linear systems and preconditioners 

Solving circulant linear system is immediate, by employing the factorization 
(12.11) . and the computation requires just two fft's. The algorithm is imple- 
mented in the functions mldivide and mrdivide, placed into the Qsmcirc 
directory, and is accessible via the usual matrix left /right division operators. 

To solve Toeplitz linear systems, one possibility is to use an iterative 
solver, either user supplied, or among those (peg, gmres, etc.) available in 
Matlab. This can be done transparently, taking advantage of the compact 
storage and fast matrix-vector product provided by the toolbox. 

It is usual to employ preconditioners to speed up the convergence of iter- 
ative methods. In the case of Toeplitz linear systems, it has been proved that 
various classes of circular preconditioners guarantee superlinear convergence 
for the conjugate gradient method; see [2]. 

The function smtcprec, included in the toolbox, provides the three best 
known circulant preconditioners, and can be easily extended to include more. 
For a given Toeplitz matrix A, the function can construct the Strang pre- 
conditioner JTE\, which suitably modifies the matrix to make it circulant, the 
so-called optimal preconditioner 0], which is the solution of the optimization 
problem 

min lie — A\\f, 
cec 

where C is the algebra of circulant matrices and || ■ ||f denotes the Frobenius 
norm, and the superoptimal preconditioner [3], [191 120] , which minimizes 

\\I-C-^A\\f 

for C G C. While the Strang preconditioner is defined only when A is 
Toeplitz, the optimal and superoptimal preconditioners can be computed 
for any matrix; the function smtcprec allows this, though the computation 
of the preconditioner is fast only for a Toeplitz matrix. 
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The code in the functions Strang, optimal, and superopt (see Table E]) 
was developed by one of the authors during the research which led to [22] , and 
the details of the algorithms are described in that paper. To compute the su- 
peroptimal preconditioner, it is possible to use either the method introduced 
in [3j, or the one from [20]. It is remarkable that, since the second algorithm 
is based on the use of certain Toeplitz matrices, our implementation is greatly 
simplified, as it performs the computation using the arithmetics provided by 
the toolbox itself. 

Using circulant preconditioners with the iterative methods available in 
Matlab is straightforward, as these functions use the matrix left division 
to apply a preconditioner, and so take advantage of the storage and fast 
algorithms furnished by our toolbox. For example, the instructions 

T=smtgallery ( 'gaussian' ,5000) ; 

b=T*oiies (5000,1) ; 

C=smtcprec ( ' Strang' ,T) ; 

[x,f lag, relres , iter] =pcg(T,b, [] , [] ,C) ; 

create a Gaussian linear system of dimension 5000 (see Section 12. 4p with 
prescribed solution, and solve it by the conjugate gradient method, precon- 
ditioned by the Strang circulant preconditioner. 



Preconditioners 


Direct solvers 


smtcprec 


circulant preconditioners 


toms729 


Toeplitz solver 


Strang 


Strang preconditioner 


tils 


Toeplitz LS solver 


optimal 


optimal preconditioner 


tsolve 


user supplied function 


superopt 


superoptimal precond. 


tsolvels 


user supplied fvmction 


General functions 


issmcirc 


true for smcirc object 


smtconf ig 


toolbox configuration 


issmtoep 


true for smtoep object 


smt const 


toolbox constants setting 


smtcheck 


check toolbox installation 


smtgallery 


test matrices 



Table 3: Computational and general functions 



Besides the iterative methods, there are also many fast and superfast 
direct solvers for a Toeplitz linear system, and some of them have been im- 
plemented in publicly available subroutines. With our toolbox, we distribute 
two of them, having computational complexity O(an^); they are called when 
one of the matrix division operators (either \ or /) are used to invert a 
Toeplitz matrix. The related files, listed in Table[3l are placed in the private 
subdirectory of smt/Osmtoep. 

The first one, toms729, is an implementation of the extended Levinson 
algorithm for nonsingular Toeplitz linear systems [lOj, written in Fortran, for 
which a Matlab MEX gateway is available [1]. This solver, which has been 
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implemented only for real matrices, calls the dsytep subroutine from [10\ if 
the system matrix is symmetric, and dgetep in the general case, with the 
pmax parameter set to 10. 

When the linear system is overdetermined (and full-rank) the toolbox 
calls the C-MEX program tils, developed in [15], which converts it into a 
Cauchy-like system Cy = f , and computes its least-squares solution as the 
Schur complement of the augmented matrix 





I 


c 





Mc = 


C* 












/ 






using the generalized Schur algorithm with partial pivoting. Complex linear 
systems are supported. If the matrix is either underdetermined or rank- 
deficient, an error is returned. 

It is possible for the user to use different algorithms, by supplying the 
functions tsolve, for a nonsingular linear system, or tsolvels, for least- 
squares, overwriting those placed in the directory smt/Osmtoep/private, 
and changing the default behaviour of the toolbox with the smtconf ig com- 
mand. For example, entering from the command line the instructions 

smtconfig intsolve off 
x=T\b 

disables the solver toms729, and solves the linear system Tx = b by the user 
supplied function tsolve. 

2.4 Other functions 

Besides the overloaded operators and methods located in the Osmcirc and 
@smtoep directories, some general functions, listed in Table [HI are placed in 
the main toolbox directory, and are directly accessible to the user. 

Among these functions, we find the two isXXX functions, which return 
logical values and check if the supplied parameter belongs to the XXX class, 
the smt check function, which verifies if the toolbox is correctly installed, and 
the function smtconst, intended to define global constants (for the moment 
only the dimension of the circulant embedding f l2.4p ). 

Of particular relevance is the smtconfig function, which modifies the 
behaviour of the toolbox for what concerns the display method for objects 
(display parameter), the use of the Toeplitz premultiplication routine, dis- 
cussed in Section 1212] (toeprem parameter), the warnings setting (warnings 
parameter), and the active Toeplitz solvers (intsolve and intsolvels pa- 
rameters, see Section [231) • For example. 
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smtconfig display compact (or off) selects compact display of objects, 
smtconf ig display full (or on) restores standard display method. 

All parameter are set by default to the on state; calling the command 
smtconfig with no parameters shows the state of all settings. 



Circulant matrices 


crrsind 


uniformly distributed random matrix 


crrcindn 


normally distributed random matrix 


Toeplitz matrices 


algdec 


matrix with algebraic decay 


expdec 


matrix with exponential decay 


gaussian 


Gaussian matrix 


tchow 


Chow matrix 


tdrsunadcih 


matrix of 0/1 with large determinant or inverse 


tgrcar 


Grcar matrix 


tkms 


Kac-Murdock-Szego matrix 


tparter 


Parter matrix 


tphcins 


rank deficient matrix from [9] 


tprand 


uniformly distributed random matrix 


tprandn 


normally distributed random matrix 


tprolate 


prolate matrix 


ttoeppd 


symmetric positive definite Toeplitz matrix 


ttoeppen 


pentadiagonal Toeplitz matrix 


ttridiag 


tridiagonal Toeplitz matrix 


ttriw 


upper triangular matrix discussed by Wilkinson 



Table 4: Test matrices available in the smtgallery function 



We end this section reporting another important feature of the toolbox. 
A collection of test matrices is available in smtgallery, which is modelled 
on Matlab gallery function, but returns structured objects. The collec- 
tion, listed in Table IH includes random matrices, three matrices studied 
in [23] (algdec, expdec and gaussian), one from [9j], and all the Toeplitz 
matrices provided by gallery, most of which come from [12]. The syntax 
of smtgallery is in the same style as the Matlab gallery function, and 
documentation is provided for each test matrix. For instance, 

T=smtgallery( 'gaussian' ,7) constructs a 7 x 7 Toeplitz Gaussian matrix, 
C=smtgallery ( ' errand ' , 7 , ' c ' ) returns a random complex circulant matrix, 
help private /tchow displays the help page for the Chow matrix. 

3 Implementation issues 

The toolbox is entirely written in the Matlab programming language. The 
version officially supported is 7.7 (i.e., release 2008b), anyway we tested it 
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on previous versions, way back to 6.5, without problems. To install it, it is 
sufficient to uncompress the archive file containing the software, creating in 
this way the directory smt and its subtree. This directory must be added to 
the Matlab search path by the command addpath, in order to be able to use 
the toolbox from any other directory. 

As noted in the previous sections, a great effort has been devoted to catch 
all possible user's errors, and to reproduce the standard behaviour of Matlab, 
for example for what concerns the output of each function in the presence of 
empty or scalar arrays in input. Since these features are scarcely documented 
in Matlab manuals, our choices are mostly due to experimental tests. 

We remark that the smtconf ig function, described in Section [2^ relies 
on the use of warnings, so issuing the Matlab command warning on restores 
the initial configuration of the toolbox, while warning off may cause un- 
predictable results. 

Some of the toolbox functions use the isfloat command, which was 
introduced in version 7 of Matlab. For those who are using version 6.5, a 
patch for this function is included in the software; see the README.txt ffie in 
the main toolbox directory. 

The two programs used to solve Toeplitz linear system are the only ones 
which need to be compiled: toms729 [10] was written in Fortran, and uses the 
Fortran-MEX gateway from pLj, while tils \L5l was originally developed as a 
C-MEX program. The MEX interface is a library, distributed with Matlab, 
which allows to dynamically link a Fortran or C subroutine to Matlab, and to 
exchange input and output parameters between the compiled program and 
the environment, using the usual Matlab syntax. 

Both Toeplitz solvers can be easily compiled under Linux, using the 
Makefile placed in the directory smt/Osmtoep/private, and we provide 
precompiled executables for 32 and 64 bits architectures. Compiling the 
same programs under Windows is a bit more involved: we used a porting 
of the GNU-C compiler [l7j and the "MEX configurator" Gnumex [16], but 
precompiled executables are available for various Matlab versions; see the 
README.txt ffie and the content of the smt/Osmtoep/private directory. 

4 Conclusion 

The Structured Matrix Toolbox is a Matlab package which implements op- 
timized storage and fast arithmetics for circulant and Toeplitz matrices, of- 
fering a robust and easily extensible framework. The toolbox is available at 
the web page http://bugs.unica.it/smt/. We are currently performing 
numerical tests to assess its performance, and there are plans to extend its 
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functionality by adding the support for other classes of structured matrices. 



References 

A. Arico and G. Rodriguez. toms729gw: a Matlah (Fortran) MEX Gate- 
way for TOMS Algorithm 729, by P. C Hansen and T. Chan. University 
of Cagliari, 2008. Available at: http : / /bugs . unica . it/~gppe/ soft/ [ 

R. H. Chan and X.-Q. Jin. An Introduction to Iterative Toeplitz Solvers. 
SIAM, Philadelphia, 2007. 

R. H. Chan, X.-Q. Jin, and M. C. Yeung. The circulant operator in the 
Banach algebra of matrices. Linear Algebra AppL, 149:41-53, 1991. 

T. F. Chan. An optimal circulant preconditioner for Toeplitz systems. 
SIAM J. Scz. Stat. Comput, 9(4):766-771, 1988. 

P. J. Davis. Circulant Matrices. Wiley, New York, 1979. 

J. Demmel and J. Dongarra. ST-HEC: Reliable and scalable software 
for linear algebra computations on high end computers. Available at: 
|http : //www . cs ■ berkeley . edu/^demmel/Sca-LAPACK-Proposal . pdf , 
2005. 

M. Frigo and S.G. Johnson. The design and implementation of FFTW3. 
IEEE Proc, 93(2):216-231, 2005. Software available at: 
http : //www . f f tw . org/[ 

P. C. Hansen. Regularization tools: A Matlab package for analysis and 
solution of discrete ill-posed problems. Numer. Algorithms, 6:1-35, 1994. 
Software available at: http://www2.imm.dtu.dk/~pch/Regutools/. 

P. C. Hansen. Rank-Deficient and Discrete Ill-Posed Problems: Numer- 
ical Aspects of Linear Inversion. SIAM Monographs on Mathematical 
Modeling and Computation. SIAM, Philadelphia, 1998. 

P. C. Hansen and T. F. Chan. Fortran subroutines for general Toeplitz 
systems. ACM Trans. Math. Software, 18(3):256-273, 1992. 

P. C. Hansen, M. Jacobsen, and T. K. Jensen. MOORe Tools: Modu- 
lar Object Oriented Regularization Tools. Technical University of Den- 
mark, Informatics and Mathematical Modelling, 2006. Available at: 



http : //www2 ■ imm . dtu . dk/~pch/MOOReTools/ 



17 



[12] N. J. Higham. The Test Matrix Toolbox for Matlab (Version 3.0). Tech- 
nical Report 276, Manchester Centre for Computational Mathematics, 
1995. 

The MathWorks, Natick. Matlab ver. 7. 7, 2008. 

J. Nagy. RestoreTools, An Object Oriented Matlab Package for Image 
Restoration. Emory University, Department of Mathematics and Com- 
puter Science, 2007. Available at: 

[http : //www .mathcs . emory . edu/~nagy/RestoreTools/| 

G. Rodriguez. Fast solution of Toephtz- and Cauchy-like least squares 
problems. SIAM J. Matrix Anal. AppL, 28(3):724-748, 2006. 

SourceForge.net. Gnumex, 2000. Available at: 
,http : //gnumex. sourcef orge .net/, 

SourceForge.net. MinGW, 2008. Available at: 



|http : / / www . mingw . org/ 



G. Strang. A proposal for Toeplitz matrix calculations. Stud. Appl. 
Math., 74:171-176, 1986. 

M. Tismenetsky. A decomposition of Toeplitz matrices and optimal 
circulant preconditioning. Linear Algebra AppL, 154/156:105-121, 1991. 

E. E. Tyrtyshnikov. Optimal and superoptimal circulant precondition- 
ers. SIAM J. Matrix Anal. Appl, 13(2):459-473, 1992. 

M. Van Barel. Software produced by members of the MaSe-Team (Ma- 
trices having Structure). Katholieke Universiteit Leuven, Department of 
Computer Science, 2008. Available at: 



http : / /www . cs . kuleuven . ac . be/~marc/ software/ 



C.V.M. van der Mee, G. Rodriguez, and S. Seatzu. Fast computation of 
two-level circulant preconditioners. Numer. Algorithms, 41(3):275-295, 
2006. 

C.V.M. van der Mee and S. Seatzu. A method for generating infinite 
positive self-adjoint test matrices and riesz bases. SIAM J. Matrix Anal. 
AppL, 26(4):1132-1149, 2005. 

The Working Group on Software (WGS). The Control and Systems 
Library SLICOT, 1998. Available at: jhttp: //www, slicot . org/] 



18 



[25] The Working Group on Software (WGS). 
Network NICONET, 1998. Available at: 



The Numerics in Control 



http : //www ■ icm . tu-bs . de/MICONET/ 



[26] The Working Group on Software (WGS). Basic Software Tools for 
Structured Matrix Decompositions and Perturbations, 2001. Available 



at: http : //www . icm . tu-bs . de/NIC0NEt/MICtasklB7html 



19 



